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Abstract The Helioseismic and Magnetic Imager onboard the Solar Dynamics Obser- 
vatory (SDO/HMI) provides continuous full-disk observations of solar oscillations. We 
develop a data-analysis pipeline based on the time-distance helioseismology method 
to measure acoustic travel times using HMI Doppler-shift observations, and infer solar 
interior properties by inverting these measurements. The pipeline is used for routine 
production of near-real-time full-disk maps of subsurface wave-speed perturbations and 
horizontal flow velocities for depths ranging from to 20 Mm, every eight hours. In 
addition, Carrington synoptic maps for the subsurface properties are made from these 
full-disk maps. The pipeline can also be used for selected target areas and time periods. 
We explain details of the pipeline organization and procedures, including processing 
of the HMI Doppler observations, measurements of the travel times, inversions, and 
constructions of the full-disk and synoptic maps. Some initial results from the pipeline, 
including full-disk flow maps, sunspot subsurface flow fields, and the interior rotation 
and meridional flow speeds, are presented. 
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1. Introduction 

The Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory (SDO/HMI:| 
ISchou et al\ 1201 11 1 observes the solar full-disk intensity, Doppler velocity, and vec- 
tor magnetic field of the photosphere with high spatial resolution and high temporal 
cadence. Similar to the Michelson Doppler Imager (MDI: [Scherrer et al.\ 1 1995b . an 
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instrument onboard the Solar and Heliospheric Observatory (SOHO), the HMI Dopp- 
lergrams are primarily used for helioseismic analysis to investigate the interior structure 
and dynamics of the Sun. Helioseismology data analysis pipelines are planned for near 
real-time analyses of the observations in order to provide the analysis results to the 
helioseismology and solar physics communities. The time-distance analysis pipeline 
is one of the pipelines for local helioseismology studies, and other pipelines include 
ring-diagram analysis and farside active region imaging. The time-distance pipeline is 
designed for routine production of nearly full-disk subsurface wave-speed perturbations 
and horizontal flow fields every eight hours, as well as synoptic flow maps for every 
Carrington rotation. It can also be used to analyze specific target areas and time periods. 

Time-distance helioseismology was first introduced by Duvall et al. (1993, 1996), 
and it has developed rapidly since then. Different inversion techniques were introduced 
and tested. The LSQR algorithm, introduced by IKosovichevI i 19961 1 and used later by 
Zha o, Kosovichev, and Duvall] ( 120011 ). solves the inversion problem in a least-squares 
sense in the spatial domain by an iterative approach. The Multi-Channel Deconvolution 
(MCD) method, introduced by Uacobsen et al. ( 1999 ) and widely used in later studies 
(e.g. ICouvidat et ail 120041 ). solves the least-squares problems in the Fourier domain. 
Later, ICouvidat et ail ([2005 ) applied a horizontal regularization procedure for this in- 
version technique. More recently, an optimally localized averaging (OLA) inversion 
scheme was introduced to study the solar subsurface flow fields (Jackiewicz, G izon, and Birch, 2008| ).l 

Different types of sensitivity kernels, which describe the relationship between the 
travel times and interior properties, were also introduced and used in the time-distance 

inversion problems . IKosovichevId 1 9961 ) first used ray -path approximation kernels, |Jensen, Jacobsen, and Christensen-D 
(12000b introduced Fresnel-zone kernels; and B irch and Kosovich ev (2000), [Birch, K osovichev, and Duvall|| 
( 120041 ). and lBirch and Gizonl (120071 ) investigated Born-approximation kernels for both 
sound-speed structures and flow fields. ICouvidat et al. I ( 120061 ) compared subsurface 
sound-speed perturbation structures inferred from these different types of kernels, and 
found that the inversion results obtained with the different kernels were basically con- 
sistent. 

Important results on the solar interior properties have been obtained from the time- 
distance studies as well as from other local helioseismology techniques. On global 
scales, poleward meridional flows were found below the photosphere (Gil es et al., 1997| ),1 
and solar-cycle dependent meridional flow variations were also investigated and dis- 
cussed ( |Chou and Dai, 2001||Beck, Gizon, and Duvall, 2002[|Zhao and Kosovichev, 2004| ).| 
On local scales, subsurface sound-speed perturbations and flow fields were derived for 
supergranulation dKosovichev and Duvall, 1997[|Duvall et a l. , 1997 ; Duvall and Gizon, 2000H 
|Sekii et al., 2007 ; Jackiewicz, Gizo n, and Birch, 2008) and for sunspots dKosovichev, Duvall, and Scherrer, 2000H 
|Gizon, Duvall, and Larsen, 2000[|Zhao, Kosovichev, and Duvall, 2001 [ICouvidat et al., 200"6[ l 
|Zhao, Kosovichev, and Sekii, 20 1 0| > . Additionally, time-distance helioseismology was 
used to detect the emergence of active regions before their appearances in the photo- 
sphere dKosovichev, Duvall, and Scherrer, 2000||Jensen et al., 2 001 ; Zha rkov and Thompson, 2008| ,[ 
to image large active regions on the farside of the Sun (Zh ao, 2007[ Ilonidis, Zhao, and Hartlep, 2009]),[ 
and to measure sound-speed perturbations in the tachocline ( |Zhao et al., 2009| ). These 
results are important for space-weather forecasting and understanding the mechanisms 
for the generation of solar magnetism. The time-distance helioseismology pipeline 
analyses, based on the high spatial-resolution and high temporal-cadence observations 
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from HMI, will greatly advance our knowledge of the interior processes and their 
connections with solar activity above the photosphere. 

However, one should keep in mind that the physics of solar oscillations in the tur- 
bulent magnetized plasma is very complicated, and that the helioseismic techniques 
are still in the process of being developed. Because of limited knowledge of the wave 
physics and complexity of the MHD turbulence, there may be systematic uncertain- 
ties in the local helioseismology inferences, particularly in strong magnetic-field re- 
gions of sunspots. For example, Lindsey and Braun (2005a, 2005b) argued that the 
outgoing and ingoing travel time asymmetries observed in sunspot areas might be 
caused by a "shower-glass effect". [Schunker et al. I (12005 ) found that the inclined mag- 
netic field in sunspot penumbrae might cause variations of measured acoustic travel 
times. IZhao and Kosovichevl (120061) found that this inclined magnetic field dependent 
effect does not exist in the measurements obtained from the MDI intensity observa- 
tions. Later, Rajaguru et al. (120061 ) found that the non-uniform acoustic power dis- 
tribution in the photosphere also contributed to measured travel-time shifts in active 
regions if a phase-speed filtering procedure was applied. This effect was then studied by 
Parchevsky, Zhao, and Kosovichev|lf2008l) and Hanasoge et al. d2008b numerically, and 



also by Nigam and Kosovichev| ( [20 10) analytically. More recently. iGizon et ai\ (12009 ) 
derived a sunspot's subsurface flow fields after applying ridge filtering, and their in- 
ferred results did not agree with the previously inverted results with a use of phase-speed 
filtering. Using high spatial-resolution observations from Hinode, Zh ao, Kosovichev, and SekH 
(2010) showed that the principal results on sunspot structure did not depend on the 
use of phase-speed filtering. However, significant systematic uncertainties in sunspot 
seismology remain and need to be understood, and these are being actively studied by 
the use of numerical simulations. For a recent review of the sunspot seismology and 
uncertainties see IKosovichevI ( 120 1 01 ) . 

Despite the ongoing discussions of accuracy of time-distance measurements and 
interpretation of inversion results, it is useful to provide the measured travel times 
and the inversion results for investigations of structures and flows below the visible 
surface of the Sun. As the flow chart of the pipeline in Figure[T]shows, we apply phase- 
speed filtering to the HMI data, compute cross-covariances of the oscillations, and use 
two different travel-time fitting procedures to derive the acoustic travel times. We then 
perform inversions using two different sets of travel-time sensitivity kernels, based on 
the ray -path and Born approximations, to infer subsurface wave-speed perturbations and 
flow fields, using the MCD inversion method. We provide online access to the measured 
travel times, full-disk subsurface wave-speed perturbations and flow maps calculated 
every eight hours. In addition, we provide synoptic flow maps for each Carrington 
rotation. In this article, we describe details of the acoustic travel time measurement 
procedure in Section [2] and the inversion procedure in Section [3] We describe the 
pipeline data products and present initial HMI results in Section [4] and summarize our 
work in Section 
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Figure 1. Flow chart for the HMI time-distance helioseismology data analysis pipeline. 



2. Acoustic Travel-Time Measurement 

2.1. Tracking and Remapping 

The SDO/HMI continuously observes the full-disk Sun, providing Doppler velocity, 
continuum intensity, line-depth, line-width, and magnetic field maps with a 45-second 
cadence, and also vector magnetic field measurements with a cadence of 12 minutes. 
Each full-disk image has 4096 x 4096 pixels with a spatial resolution of 0.504 arcsec 
pixel -1 (i.e., approximately, 0.03 heliographic degree pixel -1 at the solar disk center). 
The Doppler observations are primarily used for helioseismic studies. The observ- 
ing sequences, algorithms for deriving Doppler velocity and magnetic field, and other 
instrument calibration issues are discussed bv lSchou et q/.ld201 11 1. 

As illustrated in Figure Q] the primary input for the pipeline is Dopplergrams, al- 
though in principle, the HMI intensity grams and line-depth data can also be analyzed in 
the same manner. Users of the pipeline can select specific areas for analysis, preferably 
within 60° from the solar disk center. In practice, the users provide the Carrington 
longitude and latitude of the center of the selected area, and the mid-time of the selected 
time period, then the pipeline code selects an area of roughly 30° x 30° centered at the 
given coordinate, and for a time interval of eight hours with the given time as the middle 
point. The data for this selected area and the time period are then tracked to remove solar 
rotation, and remapped into the heliographic coordinates using Postel's projection (also 
known as azimuthal equidistant projection) relative to the given area center. Normally, 
the tracked area consists of 512 x 512 pixels with a spatial sampling of 0.06° pixel -1 ; 
and the temporal sampling is the same as the observational cadence. Cubic interpolation 
is used for the pixels not located on the observational grid. 

Figure [2] shows typical HMI k-u and time-distance diagrams obtained from one 
tracked and remapped area. The selected area covers 30° in latitude and has an apparent 
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k (Mm -1 ) 




differential rotation over this span. However, for fast computations, only one uniform 
tracking rate corresponding to the Snodgrass rotation rate at the center of this area is 
used. The Snodgrass rotation rate is: 2.851 — 0.343 sin 2 cf> — 0.474 sin 4 jurad s~\ 



where <fi is latitude (Snodgrass and Ulrich, 1990 . The uniform tracking rate of the se- 
lected area results in a differential rotation velocity in the inverted horizontal velocity 
fields. This differential rotation velocity is removed from the full-disk flow maps after 
averaging over the whole Carrington rotation, and only the residual flow fields are given 
as the final results (see Section 14.1b . But for the user-selected areas, the differential 
rotation is kept in the inversion results. 
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Table 1. Phase-speed filtering parameters used for the selected travel 
distances (annulus ranges). 



nnniiliic Nr\ 




nhasp QTiPfH 


FWHM 




(heli.ogr3.phic degree) 






l 


0.54 - 0.78 


3.40 


1 .0 


2 


0.78-1.02 


4.00 


1.0 


3 


1.08-1.32 


4.90 


1.25 


4 


1.44-1.80 


6.592 


2.149 


5 


1.92-2.40 


8.342 


1.351 


6 


2.40 - 2.88 


9.288 


1.183 


7 


3.12-3.84 


10.822 


1.895 


8 


4.08-4.80 


12.792 


2.046 


9 


5.04-6.00 


14.852 


2.075 


10 


6.24-7.68 


17.002 


2.223 


11 


7.68-9.12 


19.133 


2.039 



2.2. Computing Cross-Correlations and Fitting for Travel Times 

Each tracked and remapped Dopplergram datacube is filtered in the 3D Fourier do- 
main. Solar convection and /-mode oscillation signals are removed first, and then 
phase-speed filtering is applied following the procedures prescribed bv lCouvidat et al.\ 
(2005 ). For the travel-time measurements, for each central point we select 11 annuli 
with various radii and widths chosen from our past experience with MDI analyses. 
All of the phase-speed filtering parameters, including the central phase-speed, filter 
width, and the corresponding inner and outer annulus radii are presented in Table 1. 
The phase-speed filter is a Gaussian function of the wave's horizontal phase speed. It 
selects wave packets traveling between the central points and the annuli for the selected 
distances. This filter helps improve the signal-to-noise ratio, and also removes leak- 
age from low-degree oscillations. After the filtering, the data are transformed back to 
the space-time domain for cross-covariance computations. Two different fitting meth- 
ods are used to derive the acoustic travel times from the cross-covariances: a Gabor 
wavelet fitting ( |Kosovichev and Duvall, 1997| , and a cross-correlation method based 
on seismology algorithms ( |Zhao andJo rdan, 1998) adopted by Gizon and Birch (2002; 
GB algorithm hereafter). A detailed description of the filtering procedure, the cross- 
covariance computations, the two types of travel-time fittings, comparison of the travel 
times derived from the two fitting methods, and the measurement error estimates is 
given bv lCouvidat et al. I (1201 11 ). In particular, it has been shown there that the two defi- 
nitions of the acoustic travel times and the two fitting methods give generally consistent 
results in the quiet-Sun regions, but may give different results in active regions. The 
differences are currently not well understood, but in the pipeline we implement both 
travel-time definitions. 

In the Postel-projected maps, the exact distance between two arbitrary points cannot 
be calculated using the Cartesian coordinates. Thus, when an annulus is selected around 
a given location, some additional computations are needed to determine the exact great- 
circle distance between the points. The formula to determine the great-circle distance 
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is: 

A = arccos(sin6'i sin #2 + cos#i cos #2 cos(<^i — 02)), (1) 

where (f?i, (j>i) and (#2, ^2) are the heliospheric longitude and latitude coordinates for 
the two separate locations. 

To facilitate the inversions for subsurface flow fields, each annulus is divided into 
four quadrants representing the North, South, East, and West directions ( |Kosovichev and Duvall, 1997| .l 
So, for each annulus and each travel-time fitting method, we obtain six measurements of 
acoustic travel times, corresponding to the outgoing and ingoing, East-, West-, South-, 
and North-going directions. We then combine these travel times to obtain one map for 
the mean travel time and three maps for the travel-time differences. These travel times 
are: T mea n (average of outgoing and ingoing travel times), r j (difference of outgoing 
and ingoing travel times), r wc (difference of West- and East-going travel times), and 
r ns (difference of North- and South-going travel times). These four travel-time maps 
for each annulus are archived and available through the HMI Data Record Management 
System (DRMS). 

3. Subsurface Wave-Speed Perturbation and Flow Field Inversions 

The acoustic travel times are derived by two different fitting methods: the Gabor wavelet 
function and the GB algorithm. Then, as illustrated in Figure [T] to infer the subsurface 
wave-speed perturbations and flow velocities, the Gabor-wavelet travel times are in- 
verted using the ray-path approximation sensitivity kernels, and the GB travel times 
are inverted using the Born-approximation sensitivity kernels. The Born-approximation 
kernels are calculated based with the filter and window parameters of the GB fits. 

3.1. Inversions 

Both the ray-path and Born-approximation kernels have been used in previous time- 
distance studies (e.g. Zh ao, Kosovichev, and Duvail[l2001UCouvidat et al.\[2QQ6i . De- 
tails of the kernel calculations and their comparisons will be given in a separate paper. 

We employ the MCD inversion method ( Jacobsen et al, 1999| > with a horizontal 
regularization ( |Couvidat et al, 20051 1 . F° r the wave-speed perturbation inversions, the 
linearized equation relating the measured mean travel times and the subsurface wave- 
speed perturbations are: 

ijk 

where Ssijk is the relative wave-speed perturbation Scijk/cijk approximated by piece- 
wise constant functions on the inversion grid, and is a matrix of the discretized 
sensitivity kernel. Here, A and /i label the coordinates of the central points of the annuli 
in the observed areas, v labels the different annuli, and i, j, and k are the indices for the 
discretized three-dimensional space for inversions. Usually, the horizontal coordinates 
of the inversion grid (i and j indices) are selected at the same locations as the central 
travel-time measurement points (A and fi indices). In the first-order approximation, the 
sensitivity kernels are calculated for a spherically symmetric solar model and do not 
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depend on the position on the solar surface. Therefore, in this case Equation (O is 
actually equivalent to a convolution in the horizontal domain, which can be simplified 
as a direct multiplication in the Fourier domain: 

8t v {k\,k^) = ^ AI(k, x , K^SSk^Kx, k^), (3) 
k 

where Sf, A, and 5s are the 2D Fourier transforms of St, A, and 6s, respectively; k 
is the same as in Equation (2); n\ and k m are the wavenumbers in the Fourier domain 
corresponding to A and fi of the spatial domain. For each (k\, k^), the equation in the 
Fourier domain is a matrix multiplication: 

d = Gm, (4) 

where 

d = {(St^ka,^)}, G = {AI(k x , k^)}, m = {Ss k (n\, k m )}. 

Thus, we have a large number of linear equations describing the depth dependence of 
the Fourier components, and each linear equation can be solved in the least squares 
sense. After all these equations are solved, and m is obtained for each (k\, the 
values of dsijk are calculated by the inverse 2D Fourier transform. 

Equation © is ill-posed, and regularization is required to obtain a smooth solution. 
The regularized least-squares algorithm is formulated as: 

mm{\\(d-Gm)\\l + \ 2 ( K )\\Lm\\ 2 2 }, (5) 

where ||...||2 denotes the L2-norm, L is a regularization operator, and X(n) is a reg- 
ularization parameter. We choose L to be a diagonal matrix whose elements are the 
inverse of the square root of the spatial sampling Az at each depth. Such weighting 
is necessary because the grid in the vertical direction is chosen to be approximately 
uniform in the acoustic depth, meaning that the spatial sampling of deep layers is larger 
than the sampling of shallower layers. The regularization parameter is taken in the form 
of A 2 (k) = A 2 , + A 2 (k), where A„ and A^ are vertical and horizontal regularization 
parameters. The purpose of introducing Xh is to damp the high wavenumber compo- 
nents that may lead to noise amplification. Following the discussion of lCouvidat et al.\ 
( 120051 ). we choose A^ = A2K 2 with A2 as a constant. 

Because the regularization is applied in the Fourier domain, it is quite difficult to 
use different regularization parameters for different horizontal locations of the same 
region. Sometimes, different regularization parameters are needed because in active 
regions the noise level may be different from the quiet-Sun regions, as we discuss in 
Section 3.3. Thus, it is necessary to implement into the analysis pipeline another inver- 
sion technique, the LSQR algorithm, which solves Equation (2) in the space domain by 
an iterative approach. This implementation is currently under development. 

3.2. Inversion Depth and Validation of Inversions 

For both the wave-speed and flow-field inversions, and for both the ray-path and Born- 
approximation inversions, we select a total of 11 inversion depths as follows: 0-1, 
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Table 2. Measured mean travel times and uncertainties for both quiet Sun regions and active regions. 



annulus No. 


mean travel time 
(min) 


uncertainty for quiet regions 
(min) 


uncertainty for active regions 
(min) 


1 


11.87 


0.062 


0.17 


2 


18.82 


0.061 


0.25 


3 


21.76 


0.11 


0.26 


4 


25.85 


0.11 


0.19 


5 


28.69 


0.14 


0.15 


6 


31.18 


0.14 


0.14 


7 


35.07 


0.10 


0.10 


8 


38.86 


0.12 


0.11 


9 


42.46 


0.11 


0.093 


10 


46.63 


0.14 


0.11 


11 


50.26 


0.15 


0.11 



1 - 3, 3 - 5, 5 - 7, 7 - 10, 10 - 13, 13 - 17, 17 - 21, 21 - 26, 26 - 30, and 30 - 35 
Mm. There is a total of 1 1 depth intervals. The inversion results provide the wave-speed 
perturbations and flow velocities averaged in these layers. Due to the lack of acoustic 
wave coverage in the deep interior, the reliability of inversion results decreases with the 
depth. Thus, only inversion results for the depths shallower than 20 Mm are included in 
the pipeline output. This may change in the future when more confidence is gained in 
the deeper interior inversion results. 

In recent years, several studies have been carried out to validate the time-distance 
measurements and inversions. To validate the derived subsurface flow fields, Georgobiani et «/.[ 



(2007} and lZhao et a/. I (12007) have analyzed realistic solar convection simulations and 
found satisfactory inversion results for the shallow layers covered by the simulations. 
Validations of the wave-speed perturbation inversions based on numerical simulations 
with preset structures have also been performed. Meanwhile, numerical simulations for 
magnetic structures with flows are also underdevelopment ( |Rempel, Schiissler, and Knolker, 2009{ [ 
|Stein et ah, 201 1[ |Kitiashvili et ah, 201 1) . Validations of the time-distance helioseis- 
mology techniques will be carried out as well using these simulations. 

Cross-comparisons between different helioseismology techniques, e.g. comparing 
the mean rotation speed from our pipeline analysis with global helioseismology results, 
and comparing the subsurface flow fields with results from ring-diagram analyses, will 
also be important for the validation. 

3.3. Error Estimate 

There are two types of errors in the pipeline results: systematic errors due to our limited 
knowledge of the wave physics in the magnetized turbulent medium and simplified 
mathematical formulations, and statistical errors, which are mostly due to the stochastic 
nature of the solar oscillations ("realization noise")- Here, we only discuss the statistical 
errors. 

To estimate the errors in the inversion results, we need first to estimate the uncertain- 
ties in the measured acoustic travel times. Here, we estimate the measurement uncer- 
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Table 3. Error estimates for the relative wave-speed perturbation and horizontal velocity 
inferences in the quiet-Sun (QS) and active regions (AR). 



Depth 


wave speed for QS 


velocity for QS 


wave speed for AR 


velocity for AR 


(Mm) 




(ms- 1 ) 




(ms^ 1 ) 


0-1 


2.8 X 10" 3 


7.8 


6.3 X 10~ 3 


58.3 


1-3 


4.1 X 10" 3 


7.5 


10.9 x 10~ 3 


56.4 


3-5 


6.4 x 10" 3 


8.9 


8.7 x 10" 3 


51.1 


5-7 


4.6 x 10~ 3 


9.4 


9.7 x 10~ 3 


45.1 


7-10 


4.7 x 10~ 3 


13.1 


6.7 x 10" 3 


34.5 


10-13 


3.7 x 10™ 3 


12.9 


3.1 x 10~ 3 


28.1 



tainties following the prescriptions of Jensen et al. (2003| and lCouvidat et ai\ (12006b . 
We select 25 quiet-Sun regions, and use the rms variation of mean travel times for 
different measurement distances as an error estimate for the travel times. The estimated 
uncertainties for the Gabor wavelet fitting are given in Table [2] and the uncertainties 
obtained for the GB algorithm are similar, but slightly larger for short distances and 
slightly smaller for long distances. Active regions have different measurement uncer- 
tainties. To estimate these, we selected a relatively stable sunspot, NOAA AR 11092, 
from 2 through 5 August 2010, and we assumed that the sunspot did not change during 
this period. Then we use the rms of the travel times measured inside the sunspot as an 
error estimate. Due to the evolution of the sunspot, this approach overestimates the mea- 
surement uncertainties, but can still give us an approximate estimate of measurement 
errors. These error estimates are presented in Table 0as well. 

Inversions are then performed for the same quiet-Sun regions and the active region 
to estimate the statistical errors in inversion results. Following the same approach, the 
rms of inverted wave-speed perturbations is assumed as the statistical error. The error 
estimates for both the quiet Sun and active regions are shown in Table [3] Because 
supergranular flows are dominant in the flow fields, it is difficult to estimate errors of the 
inverted velocity for the quiet Sun by this approach. Instead, we estimate the velocity 
errors based on the rotational velocity profile, which has little change in a time scale 
of a few days. While the errors for the wave-speed perturbation inside active regions 
are roughly twice of those for the quiet-Sun regions, the velocity errors inside active 
regions sometimes are seven times larger than the errors in the quiet Sun. 



4. Data Products and Initial Results from HMI 

The time-distance data analysis pipeline is used for the routine production of nearly 
real-time full-disk (actually, nearly full-disk covering a 120° x 120° area on the solar 
disk) wave-speed perturbation and flow field maps every eight hours. These maps are 
then used to construct the corresponding synoptic maps for each Carrington rotation. 
The pipeline can also be used for specific target areas, such as active regions. In this 
section, we introduce the data products from this pipeline and some initial results from 
it. 
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Figure 3. Schematic plot showing how areas are selected for a routine calculations of the full-disk 
wave-speed and flow maps. Not all of the 25 selected areas are shown. 

4. 1 . Routine Production: Full-Disk and Synoptic Maps 

For each day of HMI observations, we select three eight-hour periods: 00:00 - 07:59 
UT, 08:00 - 15:59 UT, and 16:00 - 23:59 UT. For each analysis period, we select 25 
regions, with the central locations at 0°, ±24°, and ±48° in both longitude and latitude, 
where the longitude is relative to the central meridian at the mid-time of the selected 
period. Figure [3] shows locations of these areas on the solar disk. The total number of 
areas is 25: five rows and five columns. Due to the Postel's projection, the boundaries 
of these areas are often not parallel to the latitude or longitude lines. It is also evident 
that many areas overlap, some areas overlap twice, and some overlap four times. The 
travel times and inversion results are averaged in these overlapped areas. However, in 
the areas close to the solar limb, the foreshortening effect may become non-negligible, 
but the role of this is not yet systematically studied. Users of these maps are urged to 
be cautious when using the pipeline results in the areas close to the limb. 

For each full-disk map and each synoptic map, the East - West velocity [v x \, the 
North - South velocity [v y ], and wave-speed perturbation [5c/ c] in each depth layer are 
derived with a horizontal spatial sampling of 0.12° pixel -1 . For each of the 25 areas, the 
inversion results are first obtained in the Postel-projection coordinates, and then con- 
verted into the longitude-latitude coordinates. This coordinate conversion is basically 
the inverse procedure of transforming the observed data into the Postel-projection coor- 
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Longitude 

Figure 4. A map of horizontal flow divergence for a depth range of 1 - 3 Mm and a time period of 00:00 
- 07:59 UT 19 May 2010. The display scale is from -6.2 x 10~ 4 to 6.2 x 10~ 4 s" 1 . White areas with 
positive divergence correspond to supergranulation. Note that supergranules appear larger at high latitudes 
because of the rectangular longitude-latitude map projection. 

dinates for the travel-time measurements. Cubic spline interpolation is employed. The 
results in high latitude regions are oversampled. After the coordinate transformation, the 
overlapped areas are averaged, and the final statistical errors are estimated for the whole 
procedure starting from the travel-time measurements. This includes all potential errors 
from the interpolation and remapping. The final full-disk results are saved on a uniform 
longitude-latitude coordinate grid, so each horizontal image of the subsurface layers 
has a total of 1000 x 1000 pixels covering 120° in both longitudinal and latitudinal 
directions. This coordinate choice is convenient for merging and analyzing results, but 
unavoidably distorts maps in high latitude areas. 

Figure [4] shows an example of a full-disk map of the subsurface horizontal flow 
divergence calculated from v x and v y at the depth range of 1 - 3 Mm. The positive 
divergence areas correspond to supergranulation. Such maps at various depths with 
continuous temporal coverage can be valuable for studying the structure and evolution 
of the supergranulation. Figure |5]displays the subsurface horizontal flow fields with full 
spatial resolution for a region located at the center of the map in Figure|4] Supergranular 
flows can be easily identified. 
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Figure 5. A sample of subsurface horizontal flow fields with full spatial resolution at the depth of - 1 
Mm. This area is sampled at the center of the map shown in Figure [4] The background image shows the 
line-of-sight magnetic field measured by HMI, with red as positive and blue as negative polarities. The range 
of the magnetic field is from —80 to 80 Gauss. 



Carrington Rotation 2097 
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Figure 6. Synoptic map of large-scale horizontal flows at the depth of 1 — 3 Mm for CR 2097. The back- 
ground image is the corresponding line-of-sight magnetic field, with red as positive and blue as negative 
polarities, in the range of —50 to 50 Gauss. 
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From the full-disk wave-speed perturbation and flow maps obtained every eight 
hours, we select an area 13.2° wide in longitude, i.e. —6.6° to 6.6° from the central 
meridian, to construct the synoptic wave-speed perturbation and flow maps. Since the 
Carrington rotation rate corresponds to a shift of approximately 4.4° every eight hours, 
each location in the constructed synoptic maps is averaged roughly three times {i.e. one 
whole day). The resultant synoptic map for each depth consists of 3000 x 1000 pixels. 
Since such a map is difficult to display, we show in Figure [6] a binned-down synoptic 
large-scale flow map obtained for the depth of 1 — 3 Mm for Carrington Rotation 2097 
during the period from 20 May to 16 June, 2010. The vectors in the figure are obtained 
by averaging the flow fields in areas of 15° x 15° with a 3° sampling rate. The maps 
of this type are similar to the subsurface flow maps obtained from the ring-diagram 
analysis (Hahcr et al, 20 02). From Figure [6] it can be found that the large-scale flows 
often converge around magnetic regions. 

Normally, the full-disk flow maps and the synoptic flow maps are provided as resid- 
ual flow maps after subtracting the flows averaged over the entire Carrington rotation. 
The subtracted average maps contain the differential rotation, meridional flows, and 
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Figure 8. Subsurface flow fields of AR 1 1072 at the depth of 1 - 3 Mm, obtained from 16:00 - 23:59 UT 
23 May 2010. The background image is line-of-sight magnetic field, with red positive and blue negative. The 
image is displayed with a scale from —1000 to 1000 Gs. The arrows are displayed after a 2 X 2 rebinning, 
and the longest arrow represents a speed of 300 m s - 1 . 



possibly some systematic effects. Figure Q shows examples of the subsurface differen- 
tial rotation speed and meridional flow speed obtained by averaging the calculations 
for Carrington Rotation 2097. These results are also provided online together with the 
synoptic flow maps. 

4.2. Target Areas 

As already mentioned in Section|2] the pipeline can also be run for specific target areas 
and specific time intervals. The pipeline users are required to provide the Carrington 
coordinates of the center of the target area, and the mid-time of the time interval. 

Figure|8]shows an example of a small part (approximately 1/9) of a target area, with 
an active region, AR 11072, included in it. The subsurface flow field, at the depth of 
1 — 3 Mm, is displayed after a 2 x 2 rebinning. The displayed time period, 16:00 - 
23:59 UT on 23 May 2010, is approximately 2.5 days after the start of emergence of 
this active region that was still growing. Our results clearly show subsurface outflows 
around the leading sunspot, and some converging flows inside it. Comparing with the 
previous results of IKosovichevl ( 120091 1 and |Zhao, Kosovichev , and Seloi] (120101 1. one 
may conclude that the subsurface flow fields of active regions evolve with the evolution 
of active regions. There may be prominent outflows around sunspots during their early 
growing phase and also their decaying phase, but there may also be strong converging 
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flows during the stable period. Our continuous monitoring of the solar full-disk sub- 
surface flows may give us an opportunity to statistically study the changes of sunspot's 
subsurface flows with the sunspot evolution. 

5. Summary 

We have developed a time-distance helioseismology data analysis pipeline for SDO/HMI| 
Doppler observations. This pipeline performs acoustic travel-time measurements based 
on two different methods, and conducts inversions based on two different sensitivity 
kernels calculated in the ray -path and Born approximations. The pipeline gives nearly 
real-time routine products of full-disk wave-speed perturbations and flow field maps 
in the range of depth 0-20 Mm every eight hours, and provides the corresponding 
synoptic wave-speed perturbation and flow field maps for each Carrington rotation. 
The averaged rotation speed and meridional flow speed are also provided separately for 
each rotation. In addition to these routine production, the pipeline can also be used for 
analysis of specific target areas for specific time intervals. This data analysis pipeline 
will provide important information about the subsurface structure and dynamics on both 
local and global scales, and its continuous coverage through years to come will be useful 
for understanding the connections between solar subsurface properties and magnetic 
activity in the chromosphere and corona. With the improvement of our understanding of 
acoustic wave and magnetic field interactions, and also the measurement and inversion 
techniques, the pipeline codes will be regularly revised. 
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